/* smog_pollution_sim: This program uses regression estimates of the effect of Smog Checks
on pollution to simulate certain counterfactuals (increasing station quality, and eliminating
Smog Check entirely*/

capture log close

clear all

set more off
set matsize 2000



cd /home/work/projects/smog_check/Weather
log using smog_pollution_sim.log, replace

use pollution_star_merged

xtset county date

foreach var of varlist  firstgen* othergen* {
	qui replace `var' = (`var')/1000
}


/*Save original versions of the key variables*/

gen fpr_orig = fpr_current_firstgen_90
gen firstgen_orig = firstgen_90
gen othergen_orig = othergen_90

/* Generate day-of-week and calendar week*/

gen dow = dow(date)
gen week = wofd(date)

/*Declare the "if" conditions in a local, as well as the full set of controls*/

local conditions  if year>=1998 & inlist(county,5,12,18,22,26,53)==0 & (inlist(county,3,6,11,12,13,17,23,47,58)==0 | (year !=2002&year!=2003)) 
local controls tMin tMax prec ib1.county ib1.county#i.year  i.week


capture drop _merge

/*Now run regressions for each pollutant*/

foreach pollutant in nox co pm10{

	if "`pollutant'" == "pm10"{
	
		tempfile temp 
		save `temp'
		qui keep if `pollutant' !=.

		collapse `pollutant' tMin tMax prec (last) firstgen* othergen* fpr*, by(county year week)


		qui regress `pollutant' c.fpr_current_firstgen_90#c.firstgen_90 firstgen_90 c.fpr_current_othergen_90#c.othergen_90 othergen_90 `controls' `conditions', vce(cluster county)
		//regress `pollutant' firstgen_fprbin*_90 othergen_fprbin*_90 tMin tMax prec ib1.county ib1.county#(c.date##c.date##c.date)  i.week   `conditions', vce(cluster county)
	
	}
	else{
		qui replace `pollutant' = `pollutant'*1000
	
		qui regress `pollutant' c.fpr_current_firstgen_90#c.firstgen_90 firstgen_90 c.fpr_current_othergen_90#c.othergen_90 othergen_90 `controls'  `conditions', vce(cluster county)
		//regress `pollutant' firstgen_fprbin*_90 othergen_fprbin*_90 `controls'  `conditions', vce(cluster county)
	
	}
	
	/*Baseline prediction*/
	qui predict `pollutant'_pred if e(sample)
	
	/*Impose "ideal" STAR, targetting C-FPR*/
	qui replace fpr_current_firstgen_90 = fpr_sim_firstgen_90
	qui predict `pollutant'_sim_ideal if e(sample)
	
	
	/*Eliminate Smog Check Entirely*/
	
	qui replace fpr_current_firstgen_90 = fpr_orig
	qui replace firstgen_90 = 0
	if "`pollutant'"!= "pm10" qui replace othergen_90 = 0
	qui predict `pollutant'_sim_nosmog if e(sample)

	qui replace firstgen_90 = firstgen_orig
	qui replace othergen_90 = othergen_orig
	
	if "`pollutant'" == "pm10" {
		keep county week pm10_pred pm10_sim*
		joinby county week using `temp', unm(using)
	}


}






gcollapse nox nox_pred nox_sim* co co_pred co_sim* pm10 pm10_pred pm10_sim* fpr_orig fpr_real_sim=fpr_real_sim_90 fpr_ideal_sim = fpr_sim_firstgen_90 (sum) N re_pass directed /*if year == 1998 | year == 2002 | year == 2009*/, by(county year)

/*Save this, so we can generate the maps in a separate script*/
save fpr_mapdata, replace

decode county, gen(NAME)
drop if NAME == ""
merge m:1 NAME using county_map_db

replace year =1998 if year ==.
expand 2 if _merge == 2, gen(flag)
replace year = 2009 if flag


gen pct_chg_nox = abs(nox -nox_sim_ideal)/nox
gen pct_chg_co = abs(co-co_sim_ideal)/co

#delimit ;
spmap pct_chg_nox using county_map_coord if year == 1998, id(_ID) ndfcolor(gs13) fcolor(Reds2) clmethod(custom) eirange(0 .5)
ndlabel(No Data) legtitle(% Decline in NOx Levels) legorder(lohi) legstyle(2) clbreaks(0 .025 .05 .075 .1 .15 .2 .9) legend(position(2))
title(1998) name(nox98, replace) graphregion(color(white))
;
graph export fpr_sim_nox_map_98.png, replace width(800);

#delimit;
spmap pct_chg_nox using county_map_coord if year == 2009, id(_ID) ndfcolor(gs13) fcolor(Reds2) clmethod(custom) eirange(0 .5)
ndlabel(No Data) legtitle(% Decline in NOx Levels) legorder(lohi) legstyle(2) clbreaks(0 .025 .05 .075 .1 .15 .2 .9) legend(position(2))
title(2009) name(nox09, replace) graphregion(color(white));
graph export fpr_sim_nox_map_09.png, replace  width(800);

#delimit;
spmap pct_chg_co using county_map_coord if year == 1998, id(_ID) ndfcolor(gs13) fcolor(Reds2) clmethod(custom) eirange(0 .5)
ndlabel(No Data) legtitle(% Decline in CO Levels) legorder(lohi) legstyle(2) clbreaks(0 .025 .05 .075 .1 .15 .2 .9) legend(position(2))
title(1998) name(co98, replace) graphregion(color(white));
graph export fpr_sim_co_map_98.png, replace  width(800);


#delimit;
spmap pct_chg_co using county_map_coord if year == 2009, id(_ID) ndfcolor(gs13) fcolor(Reds2) clmethod(custom) eirange(0 .5)
ndlabel(No Data) legtitle(% Decline in CO Levels) legorder(lohi) legstyle(2) clbreaks(0 .025 .05 .075 .1 .15 .2 .9) legend(position(2))
title(2009) name(co09, replace) graphregion(color(white));
graph export fpr_sim_co_map_09.png, replace  width(800);


